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ABSTRACT 


A  hypothesis  Is  made  that  the  Galerkln  Finite  Element 
Method  (GFEM)  offers  a  viable  option  to  the  traditional 
Finite  Difference  Method  (FDM)  for  numerical  weather  predic¬ 
tion.  The  shallow  water  barotroplc  primitive  equations  are 
the  forecast  equations  for  all  experiments.  The  hypothesis 
is  tested  by  observing  simple,  analytic,  atmospheric  wave 
propagation  on  uniform  and  variable  mesh  grids.  Second,  a 
strongly  forced  solution  simulating  small  scale  nonlinear 
Interactions  is  evaluated  for  both  the  GFEM  and  FDM. 
Finally,  a  variable,  moving  grid  for  a  GFEM  model  Is 
compared  to  a  uniform,  higher  resolution  GFEM  model  for  a 
strong  vortex  In  a  mean  flow.  The  GFEM  shows  a  better 
propagation  for  simple  atmospheric  waves  and  better  predic¬ 
tion  to  a  forced  nonlinear  solution  than  the  FDM  model.  A 
moving  variable  grid  follows  an  area  of  strong  gradients 
while  not  generating  noise  in  the  transition  zone. 
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I.  INTRODUCTION 


Proper  simulation  of  atmospheric  flow  Is  the  main 
objective  of  numerical  weather  prediction.  The  foundation 
of  numerical  weather  prediction  Is  a  set  of  equations 
including  the  momentum,  thermodynamic,  moisture  and  conti¬ 
nuity  equations  for  modeling  the  atmosphere.  Through 
computer  simulation,  numerical  weather  prediction  predicts  a 
future  state  based  on  Initial  conditions  describing  the 
atmosphere.  A  successful  forecast  model  must  include  small 
scale  processes.  Transports,  conversions  and  exchanges  of 
mass,  momentum  and  energy  occurring  on  the  small  scale 
represent  important  features  which  must  be  properly  simula¬ 
ted.  Feedbacks  from  the  proper  representation  of  these 
small  features  can  sometimes  markedly  influence  the  larger 
scale  solutions.  Representation  of  the  effects  of  small 
scale  processes  can  be  accomplished  directly  through 
increased  spatial  resolution.  Continued  increases  in  the 
spatial  resolution  will  eventually  allow  the  desired  process 
to  be  resolved  properly.  However,  a  doubling  of  resolution 
generally  requires  an  eight  fold  increase  in  computational 
effort.  The  value  of  Increased  spatial  resolution  must 
ultimately  be  measured  by  Its  contribution  to  the  overall 
forecast.  A  level  of  confidence  In  the  forecast  must  be 
achieved  that  the  process  Is  resolved  near  that  particular 
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grid  resolution.  The  grid  resolution  could  be  uniformly 
fine  In  the  forecast  domain.  The  computational  effort  for 
uniform  fine  mesh  models  could  be  far  In  excess  of  available 
resources.  An  alternative  to  a  uniform  fine  mesh  Is  a 
variable  mesh  where  the  fine  mesh  covers  only  regions  of 
Interest  or  high  activity. 

The  conventional  numerical  weather  prediction  forecast 
scheme,  the  finite  difference  method,  approximates  the 
partial  differential  equations  with  a  truncated  Taylor 
series.  It  has  performed  admirably  when  forecasting  the 
larger  scales  of  motion.  Technological  Improvements  In 
computing  power  coupled  with  better  understanding  of  the 
atmosphere  now  allow  the  smaller  scales  to  be  forecast. 

An  alternate  method  for  numerical  weather  prediction, 
the  GFEM  approximates  the  partial  differential  equations 
while  minimizing  the  error  between  the  actual  equations  and 
their  approximation.  This  best  fit  logically  leads  one  to 
the  expectation  that  the  GFEM  will  better  model  the  smaller 
scales  than  the  finite  difference  schemes.  This  research 
will  demonstrate  practical  aspects  of  the  GFEM  theoreti cal ly 
possible. 

In  this  dissertation,  the  GFEM  will  be  evaluated  to 
determine  Its  potential  to  model  atmospheric  flow.  Equiva¬ 
lent  GFEM  and  FDM  models  will  be  utilized  to  compare  the  two 
methods.  Simulations  of  small  scale  processes  explicitly 
resolved  on  uniform  and  variable  grids  will  be  Investigated. 


Rigorous  experiments  will  demonstrate  both  GFEM  and  FDH 
responses  for  forcing  near  the  grid  length  scale.  Finally* 
a  demonstration  will  be  made  of  a  GFEM  moving  variable  mesh 
which  moves  with  the  small  scale  process  or  feature,  and 
thereby  allows  the  fine  mesh  to  resolve  the  highly  active 
region  Initially  and  throughout  the  forecast  period. 


II.  HYPOTHESIS 


Numerical  weather  prediction  has  steadily  Improved  In 
the  simulation  of  atmospheric  flow.  Large  scale  flows  have 
been  adequately  represented  by  finite  difference  models  for 
several  decades.  Gal erkl n-type  formulations  (Cullen,  1974b; 
Hlnsman,  1975;  Stanlforth  and  Mitchell,  1978;  Cullen  and 
Hall,  1979;  Stanlforth  and  Daley,  1979;  MacPherson  and 
Aksel,  1980;  and  Sasaki  and  Reddy,  1980)  have  been  shown  to 
be  competitive  with  finite  difference  models,  but  have  not 
shown  a  marked  Improvement.  Comparing  the  current  opera¬ 
tional  models  at  selected  large  computer  centers,  one  finds 
that  two  are  Galerkln  and  three  are  finite  difference 
(Galerkln:  National  Meteorological  Center  and  Canadian 
Meteorological  Center;  and  finite  difference:  Fleet 
Numerical  Oceanography  Center,  Air  Force  Global  Weather 
Center  and  European  Center  for  Medium  Range  Weather 
Forecasting.  However,  all  centers  have  ongoing  research 
with  Galerkln  models  and  there  are  Indications  that  these 
will  give  better  long  range  forecasts.  The  dichotomy  arises 
because  the  Galerkin  applications  have  not  vindicated  them¬ 
selves  with  a  marked  Increase  In  accuracy,  but  rather  have 
shown  equivalent  accuracies.  Stanlforth  and  Mltchel  (1977) 
stated  that  ultimately  the  best  global/hemispheric  models 
will  be  a  spectral  model.  The  spectral  model  Is  based  on 
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the  Galerkln  procedure  end  the  use  of  trigonometric  func¬ 
tions  Is  particularly  appealing  for  hemispheric  or  global 
grids.  However,  where  non-uniform  grids  are  required,  the 
GFEM  Is  a  more  logical  choice. 

Narked  Improvements  In  regional  forecasts  will  not  come 
from  better  large-scale  models,  but  rather  from  models  that 
can  simulate  smaller  atmospheric  features.  The  National 
Heather  Service  Limited  Fine  Mesh  Model  is  an  improvement 
over  the  hemispheric  model,  partly  because  it  has  higher 
resolution  and  can  resolve  smaller  phenomena.  These  smaller 
features  can  also  affect  the  large  scale  flow.  Numerical 
meteorologists  have  realized  this  for  years  and  have 
attempted  to  model  these  features. 

The  GFEM  has  the  potential  to  increase  efficiently  the 
spatial  resolution  for  the  purpose  of  simulating  accurately 
the  small-scale  processes.  If  a  variable  mesh  is  to  be 
employed,  then  it  should  be  evaluated  by  proper  simulation 
of  a  small-scale  feature.  In  previous  research  (Staniforth 
and  Mitchell,  1978),  the  refined  grid  was  tested  by  pro¬ 
pagating  synoptic-scale  waves  into  the  finer  grid  and 
demonstrating  that  the  wave  could  move  into  the  finer  grid 
without  generation  of  significant  noise.  These  conditions 
represent  a  prerequisite.  If  synoptic-scale  waves  cannot 
move  freely  Into  and  through  the  variable  grid,  then  advec- 
tlon  Interactions  will  not  occur  properly  In  the  fine  grid. 
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However,  one  must  also  show  that  a  smaller  scale  atmospheric 
feature  can  be  properly  resolved  or  developed  In  the  fine 
mesh  area.  This  will  be  a  milestone*  The  efforts  of 
Stanlforth  and  Mitchell  (1978)  and  Stanlforth  and  Daley 
(1979)  have  stressed  the  movement  of  synoptic  scale  features 
Into  the  fine  mesh  area  but  have  not  demonstrated  Improved 
resolution  of  smaller  scale  atmospheric  features  In  the  fine 
mesh  area. 

The  hypothesis  Is  that  the  GFEM  Is  a  viable  option  for 
numerical  weather  prediction  when  simulating  atmospheric 
flow  on  variable  grids.  Three  separate  features  of  the  GFEM 
will  be  explored.  Each  feature  will  establish  the  creden¬ 
tials  of  the  GFEM  as  a  viable  option.  Each  feature  Is 
intimately  related  to  the  proper  representatl on  of  a  small- 
scale  phenomenon.  The  cost  effectiveness  of  the  GFEM,  when 
evaluated  In  these  contexts,  will  provide  a  measure  of  the 
potential  contribution  of  the  method. 

A.  FEATURES 

The  following  three  specific  features  will  be  explored 
to  support  the  hypothesl s. 

1 .  Variable  Grid 

Investigations  of  a  suitable  alternative  to  the 
finite  difference  models  for  a  variable  grid  will  be 
performed.  Two  basic  subdivisions  are  ava11ab1e--tr1angular 
or  rectangular.  Stanlforth  and  Mitchell  (1978)  employed  the 


variable  rectangular  grid  as  shown  In  Figure  1.  The  grid 
contained  some  areas  of  finer  resolution  which  were  not  in 
the  verification  area.  Two  points  associated  with  having 
variable  resolution  In  peripheral  regions  arise: 

(a)  Unnecessary  computational  overhead  Is  required; 
and 

(b)  Undesirable  phase  changes  occur  as  the  wave 
propagates  In  the  peripheral  regions. 

Older  (1981)  and  Woodward  (1981)  developed  a  transformation 
procedure  to  vary  smoothly  the  resolution  for  a  channel 
domain  from  a  coarse  to  a  fine  area  In  a  triangular  subdivi¬ 
sion.  A  uniform  equilateral  subdivision  is  shown  in  Figure 
2.  The  use  of  triangles  allows  the  increased  resolution  to 
occur  only  where  desired  and  not  in  peripheral  regions. 

Two  possible  choices  of  triangles  include  the  right  and  the 
equilateral.  Hinsman  (1975)  utilized  equilateral  triangles 
while  Cullen  (1974b)  utilized  near-equilateral  triangles.  Both 
reported  excellent  wave  propagation.  Kelley  and  Williams 
(1976)  utilized  right  triangles  and  experienced  very  noisy 
solutions.  Woodward  (1981)  duplicated  Kelley  and  Williams 
effort  with  equilateral  triangles  and  found  a  major 
reduction  in  the  noise. 

The  differing  geometries  and  results  thus  far  reported 
mandate  a  further  review  of  triangular  subdivisions.  It  Is 
not  obvious  which  subdivision  is  most  suitable.  A  distinct 
advantage  of  the  rectangular  subdivision  Is  that  It  allows 
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Fig.  1.  Rectangular  subdivision  of  the  Northern  Hemisphere 
on  a  polar  stereographic  projection  (Stanlforth 
and  Mitchell,  1978). 
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Fig.  2.  Triangular  subdivision  for  a  channel  In  Cartesian 
coordinates  (Woodward,  1981). 


algorithms  to  be  developed  which  take  full  advantage  of 
vector  processors.  However,  both  subdivisions  afford  the 
luxury  of  obtaining  a  variable  grid.  Therefore,  similar 
experiments  performed  on  each  type  subdivision  will  point 
out  the  advantages  and  disadvantages  of  each. 

2.  Numerical  Simulation  of  Physical  Processes 

The  variable  grid  chosen  as  the  most  suitable 
alternative  will  be  tested  while  simulating  a  physical 
process  to  show  that  a  small-scale  atmospheric  feature  can 
be  properly  portrayed  in  the  fine  mesh  area.  As  stated  in 
Chapter  II,  the  ability  of  variable  grids  to  resolve  small 
scale  atmospheric  features  has  not  been  rigorously  tested. 
The  ability  to  move  synoptic  scale  features  into  a  finer 
mesh  is  a  prerequisite  and  must  be  shown.  The  main  crux  of 
the  problem  is  what  is  happening  in  the  fine  mesh  area.  A 
proper  scheme  must  resolve  a  small  scale  feature  near  the 
smallest  grid  length.  Schoenstadt  (1980)  and  Williams 
(1981)  indicated  that  the  most  responsive  schemes  were  eithe 
a  staggered  finite  element  grid  with  primitive  variables  or 
an  unstaggered  finite  element  grid  with  vorti ci ty/di vergence 
formulation.  This  research  will  utilize  the  latter. 

The  simulated  process  will  be  that  of  a  mass  source 
analogous  to  that  found  In  the  upper  atmosphere  above  a 
hurricane  or  on  the  leading  side  of  a  strong  trough.  The 
source  will  appear  In  the  continuity  equation.  A  known  wave 


form  will  be  assumed  and  the  required  source  to  support  the 
wave  form  will  be  analytically  derived.  The  expression  for 
the  wave  form  will  allow  control  of  the  scale  of  the  process 
so  that  a  near-grid  scale  phenomena  can  be  simulated. 

3.  Moving  Grid 

The  ability  to  move  the  fine  grid  so  that  it 
remains  centered  on  an  atmospheric  circulation  will  also  be 
demonstrated.  This  capability  will  further  enhance  the 
applicability  of  the  GFEM  for  atmospheric  simulation.  If  it 
is  shown  that  small-scale  forcings  are  properly  reflected  in 


the  flow,  and  the  grid  can  be  moved  with  the  forcing 
disturbance,  then  the  GFEM  has  great  potential  for 
atmospheric  prediction. 


III.  THE  DEVELOPMENT  OF  THE  GFEM 


The  hypothesis  in  Chapter  II  proclaims  the  GFEM  as  a 
viable  option  for  simulating  atmospheric  flow.  An  under¬ 
standing  of  the  history  and  principles  of  the  GFEM  will  give 
further  insight  for  the  expected  improved  performance  as 
compared  to  the  conventional  FOM. 

The  origin  of  the  GFEM  can  be  traced  to  the  seventeenth 
century  work  by  Leonhard  Euler.  His  work  established  the 
branch  of  mathematics  known  as  calculus  of  variation.  He 
recognized  that  there  existed  a  partial  differential  equa¬ 
tion  (POE)  associated  with  the  minimization  of  a  functional. 
This  POE  has  been  appropri ately  named  the  Eul er-Lagrange 
equation.  Solving  the  PDE  was  equivalent  to  determining  a 
stationary  value  of  the  functional.  In  the  late  nineteenth 
century.  Lord  Rayleigh  furthered  the  variational  calculus  by 
representing  the  dependent  variable  as  a  mathematical 
expression,  typically  as  a  power  series.  The  stationarity 
of  the  solution  allowed  determination  of  the  unknown  coeffi¬ 
cients  of  the  power  series.  His  work  was  generalized  in  the 
early  twentieth  century  by  W.  Ritz  and  the  procedure  has 
since  been  referred  to  as  the  Rayl ei gh-Ri tz  method.  A 
limitation  in  the  solution  by  the  Rayl ei gh-Ri tz  method  is 
the  fullness  of  the  matrix  to  be  inverted.  Shortly  after 
Ritz  completed  his  work,  Galerkin  developed  a  procedure 


25 


il 

i 

i  ■ 

r 

]  using  weighted  residuals  which  had  a  wider  application  than 

classical  variational  calculus.  The  residual  Is  made  ortho- 

i  gonal  to  a  function,  called  the  test  function,  and  thereby  is 

minimized.  While  the  Rayl ei gh-Ri tz  procedure  applies  to  the 
functional  of  an  associated  Eul er-Lagrange  equation,  the 
Galerkin  procedure  pertains  to  any  POE.  The  methods  produce 
identical  results  when  applied  to  equivalent  extremum 
probl ems . 

{  The  Galerkin  procedure  was  unknowingly  well  established 

i  by  the  end  of  World  War  II.  The  rapid  technological 

advances  made  during  and  after  the  war,  coupled  with  the 
emergence  of  the  computer,  led  the  aircraft  industry  to 
develop  new  numerical  procedures  for  solving  stress  problems 
in  aircraft  design.  A  successful  technique  was  developed 
and  mathematicians  realized,  long  after  the  fact,  that  the 
Galerkin  procedure  was  being  utilized. 

•  The  GFEM  is  a  commonly  used  subset  of  the  set  of 


Galerkin  procedures.  After  subdivision  into  a  set  of  ele¬ 
ments,  the  domain  resembles  a  completed  jigsaw  puzzle,  and 
hence  leads  to  the  terminology  “finite  elements.”  Figures  1 
and  2  are  samples  of  such  subdivisions.  The  dependent 
variables  of  the  PDE  are  represented  as  a  linear  combination 
of  known  functions,  which  are  usually  low  order  polynomials. 
The  same  function  is  employed  as  the  test  function.  When 
substituted  into  the  PDE,  these  approxi mati ons  leave  a 
residual.  Minimizing  the  residual  completes  the  procedure. 
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The  functions  are  globally  zero  except  where  the  dependent 
variables  are  defined  near  each  individual  node.  The  pro¬ 
duct  of  functions  remains  zero  except  where  the  representing 
and  test  functions  are  both  non-zero.  This  makes  the  method 
attractive  for  computer  implementation.  It  removes  the 
matrix  “fullness"  problem  found  with  the  Rayl ei gh-Ri tz 
approach.  The  minimization  process  lies  at  the  heart  of  the 
expected  improved  performance.  Each  term  of  the  PDE  has 
been  simultaneously  approximated  and  the  error  in  those 
approximations  has  been  minimized. 
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IV.  MODEL  DESCRIPTION 

Two  different  barotropic  shal  1  ow- water  models  will  be 
employed  to  test  the  hypothesis  In  Chapter  II.  Each  model 
will  have  the  Identical  domain,  boundary  and  Initial 
conditions.  The  difference  will  lay  with  the  partial 
differential  equation  approximation.  In  one  model,  a 
Galerkin  approximation  to  the  partial  differential  equation 
will  be  used,  while  the  other  model  will  use  a  finite 
difference  approximation.  Two  versions  of  the  Galerkin 
model  will  be  evaluated.  The  first  version  will  have  trian¬ 
gular  subdivisions  and  basis  functions,  while  the  other  will 
have  rectangular  subdivisions  and  basis  functions. 

A.  GALERKIN  FINITE  ELEMENT  MODEL 

The  system  of  equations  referred  to  as  the  shallow- 
water  equations  consists  of  three  equations  with  three 
forecast  variables  u  and  v.  The  equations  are  written 
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Here  4  Is  the  geopotential  height, 

u  is  the  east/west  component  of  the  wind, 
v  is  the  north/south  component  of  the  wind,  and 
f  is  the  Coriolis  parameter 

By  expanding  4  into  a  mean  («)  and  a  deviation  (<*')  the 
equations  can  be  written 

ff  +  *D  +  +  ~(v*')  «  0  (4-4) 

il+it'  +  i£.vq,0  (4-5) 

at  ax  ax  w  v  ' 

|X  +  il’  +  H+uq*o  (4-6) 

at  ay  ay  y  ol 

here  D  is  the  divergence 

K  is  the  kinetic  energy  (per  unit  mass),  and 

Q  is  the  absolute  vorticity. 

The  primes  will  be  dropped  for  the  rest  of  the  paper  for 
cl ari ty . 

Cullen  and  Hall  (1979)  showed  that  the  accuracy  of  the 
GFEM  solution  was  better  for  the  vort i c i ty-di vergence 
formulation  of  the  shallow-water  equations  than  for  an 
increase  in  resolution  with  the  primitive  formulation. 
Williams  and  Schoenstadt  (1980)  noted  that  staggered 
variable  formulation  of  the  primitive  equations  and  the 
unstaggered  vortl cl ty-di vergence  formulation  gave  the  best 
treatment  of  geostrophic  adjustment  for  small-scale 


The  vortlcl ty-dl vergence  for*  of  the  shal low-water 
equations  also  allows  the  use  of  a  seal-implicit  tine 
scheme.  This  scheme  artificially  slows  the  propagation 
speed  of  the  fastest  gravity  waves,  which  allows  a  much 
larger  time  step  than  one  could  expect  for  a  normal  Courant- 
Fredrl ch-lewy  (CFL)  stability  criterion.  This  scheme  thus 
offsets  some  of  the  extra  computational  expense  required  to 
solve  the  system  of  equations  assembled  at  each  time  step. 

The  vortlcl ty/dl vergence  form  of  the  equations  Is 

§|  +  *0  ♦  jVlu*)  +  ay(vf)  -  0  (4-7) 


If  +  *  0  (4'8> 

j|  *  -  •^■(*0)  +  ^uq)  *  0  (4-9) 

Here  Is  the  laplaclan  operator  and  c  Is  the  -  jy) 
relative  vorticlty. 

The  velocity  can  be  written  as  the  sum  of  the  rotational  and 
Irrotatlonal  components  as 

V  ■  V.  +  V 

'4»  -X 

where  V  ■  KXV\|»  and  V  »  ?x 

r  X 
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resulting  In 
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H»2*  -  -  ^r(uQ)  -  ^r(vQ) 
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|-rV2x  ♦  v2$  ■  — (vO)  -  -i-(uO)  3  5 K  a  aK 
3t  *  ?  3xVVWi  ay'-”'"  ax  a7  ay  ay 

Equation  (4-12)  can  be  further  manipulated 

*  ♦>  ■  -  If)  -  i7«»Q>  *  fy> 


(4-10) 


(4-11) 


(4-12) 


(4-13) 


The  domain  of  Integration  is  a  channel  with  east-west 
periodicity.  The  boundary  condition  at  a  wall  is 
V  .  N  -  0 

where  N  Is  an  outward  pointing  normal  vector.  Along  the 
northern  and  Southern  walls,  the  v  component  Is  equal  to 
zero,  so  that  the  v  equation  of  motion  (4-3)  reduces  to 
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The  zonal  and  meridional  components  of  the  wind  can  be 
written  as  non-dl vergent  and  Irrotational  components  as 


and 


u  .  .  it  +  i* 
ay  ax 


v  .  +  il 

ax  ay 


Then,  along  the  north/south  walls  where  v  equals  zero,  the 
boundary  condition  Is 

♦  lx  -  o 

3X  3y  (4-14) 
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The  above  condition  Is  Imposed  by  setting 
*  «  a  constant 

when  solving  the  vorticity  equation  and 


(4-15) 


iX.  »  o  (4-16) 

ay 

when  solving  the  divergence  equation.  This  Is  an  overspe¬ 
cification  but  (4-14)  would  be  difficult  to  apply.  The 
Initial  conditions  which  will  be  presented  later  are 
specifically  selected  to  satisfy  (4-15)  and  (4-16). 

The  seml-impl 1 c 1 t  scheme  Is  Implemented  by  evaluating 
all  the  terms  on  the  left  hand  side  of  the  equations  as  an 
average  at  time  levels  (t  +  At)  and  (t  -  At)  or  with  a 
centered  time  difference  as  appropriate.  All  the  terms  on 
the  right  hand  side  are  evaluated  at  time  level  t.  The 
equations  become 

||+  *(v2x(t+At)  +  v2x(t-At))  *  -  ~-(uQ )  -  ^-(vQ)  (4-17) 

72  &  *  -  -  £(vQ)  (4-18) 

+  4(t+At)  +  4(t-At>)  -  (vQ)  -  f£)  -  ■— ( (uQ)  +  f£) 

(4-19) 


The  divergence  equation  (4-19)  can  be  solved  for  x  at 
(t  +  At)and  substituted  Into  the  equation  (4-17)  to  yield 
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-  *  n  +  ^Uv)) 


(4-20) 
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where  the  overbar  denotes  an  average  of  t  +  At  and  t  -  At. 

The  system  of  three  equations  In  three  unknowns  has  been 
reduced  to  two  Poisson  equations  and  one  Helmholtz  equation 
to  be  solved  at  each  time  step. 

The  solution  procedure  Involves  solving  the  $  equation 
(4-20)  for  a  new  i.  The  divergence  equation  (4-19)  Is  then 
solved  for  5  +  and  by  subtraction  for  |^.  Finally,  the 
vortlclty  equation  (4-18)  Is  solved  for  There  are  two 

options  as  to  which  variables  will  be  history-carrying: 
either  <j>,  ♦  and  x  or  4,  u  and  v.  The  choice  was  made  to  use 
u  and  v  as  history-carrying  variables.  They  are  updated 
after  each  time  step  following  Staniforth  and  Mitchell 
(1977)  by 

<fr(t+At)  »  2$  -  $ ( t- A t )  (4-21) 

u ( t+A t )  *  2At(-^  ft  ■  37  +  (4-22) 

*<*♦“>  ■  ft  +  A  ft1  *  v(t-at)  (4-23) 

Implementation  of  the  GFEM  Is  accomplished  as  described 
in  Chapter  III.  The  north  and  south  latitudes  are  input 
parameters,  and  the  north/south  distance  is  subdivided  Into 
N  equal  parts,  where  N  Is  also  an  Input  parameter.  The 
east/west  Increment  Is  a  function  of  the  north/south  incre¬ 
ment,  as  explained  in  the  section  describing  the  individual 
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experiments.  The  Coriolis  parameter  Mill  be  a  constant 
equal  to  the  mid-channel  value. 

The  triangular  and  rectangular  domain  subdivisions  Mill 
be  demonstrated  In  the  experiment  sections.  An  appropriate 
approximating  function  for  the  rectangular  subdivision  Is  a 
bilinear  function.  In  either  case,  the  forecast  equations 
In  Galerkln  form  are 

J  •  - «  W)V1 


*/<l7((uqVj  *  wVj))v  1  +  + 


t/n,S(t-ltlv4vi 


(4-2<) 


/»*  H  jVi  *  -  -  R(<,» jVVl 


(4-25) 


/,2(A  *  Mfwij'j  -  If  jv)vi 


‘  /lay(tu01jvj 


4  If  JvJ))vt 


(4-26) 


Here  the  Integral  sign  Implies  an  area  integral  over  the 
domain,  the  j  subscript  denotes  Einstein  summation  for  the 
dependent  variables  and  the  i  subscript  Is  the  ith  nodal 
equation . 
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The  strength  of  the  6FEM  Is  that  the  spatial  deriva¬ 
tives  of  the  dependent  variables  become  derivatives  of  the 
known  approximating  function  and,  therefore,  are  known 
exactly  at  each  node.  Vihen  using  piecewise  continuous 
linear  functions,  the  first  derivatives  are  piecewise 
discontinuous  and  second  derivatives  are  not  defined. 
Therefore,  the  second  derivatives  must  be  handled  in  a 
different  fashion, 

/’Vi)- 


/ 


57Vi>  -/«J 


3X 


;  4)  — V  V  i*  — i 

T  j¥1  J  ax'3  3x  ' 


)  ayvayJ  1J  J  dyJ 


1*  i 
ay 


— v . 
ay  i 


where  the  J>  Implies  a  line  Integral  aldng  the  east/west 
or  north/south  boundaries.  The  east/west  line  Integrals 
are  zero  because  of  perl odi ci ty, but  the  north/south  line 
Integrals  add  additional  terms  to  the  equations  and  are 
satisfied  by  the  boundary  conditions.  The  Laplaclan  terms 
become 


35 


/»*Vi  *  -Mv i  h*i  -/^  £vi  +Mvjvi 


This  Integration-by-parts  procedure  allows  the  use  of  a 
linear  function  while  approximating  second  derivatives.  The 
final  form  of  the  forecast  equations  is 
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(4-28) 
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Here  the  Helmholtz  equation  (4-27)  has  utilized 


V2X  *  DIVERGENCE  *  +  & 

and 

ii  *  -  u  f 
ay  o 

along  the  wal 1 s. 

The  line  Integral  along  the  north  and  south  boundary 
has  been  dropped  from  the  vorticity  equation  (4-28),  since 
the  value  of  4*  on  the  boundaries  is  a  constant  and  therefore 
the  value  of  |-jr  is  zero.  The  line  integral  along  the  north 
and  south  boundary  has  been  dropped  from  the  divergence 
equation  (4-29),  because  the  value  of  the  normal  derivative 
along  the  north/south  walls  is  zero,  as  stipulated  by  the 
boundary  conditions.  The  initial  conditions  will  also 
satisfy  the  condition  that  the  normal  derivative  of  x  is 
zero  along  the  north/south  walls. 

B.  FINITE  DIFFERENCE  MODEL 

The  comparison  model  to  the  GFEM  model  is  an  adaptation 
of  the  staggered,  primitive-equation  model  as  described  in 
Section  7-4  of  Haltiner  and  Williams  (1980).  The  equations 
are  in  the  flux  form  and  employ  Scheme  C  as  shown  in  Section 
7-3  of  Haltiner  and  Williams  (1980).  The  baroclinic  model 
described  there  has  been  coded  to  p? -form  In  a  barotropic 
mode. 
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Since  the  GFEM  Is  written  in  differentiated  form, 
vortlclty  and  divergence  have  been  diagnosed  from  the  primi¬ 
tive  variables  in  the  finite  difference  model  to  allow 
appropriate  comparisons. 

C.  NUMERICAL  METHODS 

Although  no  attempt  Is  made  to  optimize  computational 
efficiency  in  this  research  model,  the  GFEM  must  utilize 
efficient  numerical  techniques  to  be  considered  a  viable 
option  for  numerical  weather  prediction.  Various  solution 
procedures  are  available  for  the  GFEM  system  of  equations. 
Successive  over-relaxation  was  employed  during  the  early 
research  stages.  Several  problems  arose  which  indicated  a 
need  for  a  better  solution  procedure.  The  successive  over¬ 
relaxation  procedure  resulted  in  a  bias  if  the  sweeping 
during  each  pass  was  in  the  same  direction.  Alternating  the 
direction  alleviated  that  particular  problem,  but  the  number 
of  passes  required  to  achieve  a  desired  level  of  accuracy 
became  prohibitive.  A  second  approach  employed  a  direct 
solver  using  a  Gaussian  elimination  procedure.  The  matrices 
from  the  Galerkin  procedure  are  decomposed  into  upper  and 
lower  block  tri-diagonal  matrices.  A  preprocessi ng,  repre¬ 
senting  the  forward  substitution  stage,  can  be  done  once. 

The  back  substitution  must  be  done  each  time  a  solution  is 
desired.  The  coefficients  necessary  for  this  back  substitu¬ 
tion  are  stored  in  an  efficient  manner.  The  particular 
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algorithm  is  called  a  "skyline"  solver  and  is  described  in 
detail  in  Bathe  and  Wilson  (1976).  Skyline  refers  to  the 
compact  method  of  storing  only  those  coefficients  required. 
This  method  has  thedesired  level  of  accuracy  and  a  high 
degree  of  computational  efficiency.  Staniforth  and  Mitchell 
(1977)  employed  another  direct  solver  method  called,  the 
conjugate-  gradient  method.  This  method  is  very  computa¬ 
tionally  efficient,  even  when  the  non-constant  coefficient 
Helmholtz  equation  is  involved.  The  theoretical  operation 
count  indicates  that  for  very  large  domains  with  many 
degrees  of  freedom,  the  conjugate-gradient  method  will  be 
much  better  than  the  other  direct  solver  mentioned  above. 

The  Gauss  elimination  procedure  was  utilized  instead  of  the 
conjugate-gradient  method  for  reasons  of  expediency. 

Numerical  integration  of  the  three  forecast  equations 
involves  solving  first  a  Helmholtz  equation  for  $,  second  a 
Poisson  problem  for  ^  and  finally  a  Poisson  problem  for  x- 
The  boundary  conditions  and  the  equations  are  well  posed  for 
the  first  two  equations.  However,  the  normal  derivative  of 
X  is  equal  to  zero  for  the  last  Poisson  equation  and 
requires  special  attention.  The  solution  plus  a  constant  is 
also  a  solution.  Since  only  derivatives  of  x  are  required, 
the  shape  of  the  field  is  much  more  important  than  the 
actual  value  of  x*  To  avoid  computer  round  off  errors,  the 
average  value  of  the  terms  on  the  right  hand  side  of  (4-29) 
is  removed  at  each  time  step.  This  sets  the  constant  to 


zero  and  does  not  allow  the  solution  to  grow  by  a  constant 
each  iteration. 

The  current  technology  of  large  mainframe  computers 
incorporates  vector  processing.  Utilization  of  this 
ability  can  be  achieved  when  the  system  of  equations  can  be 
vectorized.  The  GFEM  is  derived  in  vector  formulation  and 
is  easily  vectorized. 

0.  INITIAL  CONDITIONS 

1 .  Simple  Atmospheric  Waves 

The  initial  conditions  must  allow  a  relative 
amount  of  control  for  the  desired  input  parameters  as  well 
as  satisfy  the  boundary  conditions.  The  forecast  model 
history-carrying  variables  are  4>,  u  and  v.  The  analytic 
expression  for  the  streamf uncti on ,  ip ,  is 

*  =  A  sin2  if  sin  -  *j(y-ymid)  +  f-  (A-3o) 

o 

where  A  -  amplitude  of  the  perturbation, 

W  =  width  of  the  channel, 

L  *  length  of  the  channel, 
n  *  wave  number, 

0  =  mean  flow  speed,  and 
ymid  =  middle  point  of  the  channel 
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The  first  term  in  the  expression  is  the  perturbation. 
The  second  term  represents  the  north/south  slope  necessary 
to  support  a  mean  flow  of  0.  The  third  term  is  the  mean 
depth  term.  The  geopotential  height,  $,  is  related  geostro 
phically  to  the  streamfuncti on ,  ^ ,  by 


4>  a  f0tl> 


(4-31) 


By  use  of  a  tri gonometri c  identity,  the  streamfuncti on  is 
written  as 


♦  ’  I"  -  cos^jsln^x  .  jj(y.,  ,  *  f 

0 

The  vorticity  is  given  by 

?=  2a!2A  sina2x  cos  2ojy  -  a22^-sina2> 

+  a 22  ^.Sina2X  COS  2ajy 


where 


± 


a2  = 


2irn 

L 


(4-32) 


(4-33) 


The  divergence  is  determined  through  a  linearized  form  of 
the  quasi -geostrophi  c  divergence  equation  from  Chapter  3-2 
of  Haltiner  and  Williams  (1980),  in  the  form 


V20 


-a  0^-72, 

9  3X  ^ 


(4-34) 
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The  right-hand  side  of  this  equation  becomes 


Assuming  a  solution  form  for  divergence  on  the  left-hand 
side 

0  *  cos  <»2x(Cicos  2cqy  +  C2) 
the  constants  Cj  and  C2  become 

f  f  2 

Ci  *  -  (-|-U2a.12Aa2  +  -fUa23f')  /  ( ^oq2  +  a22  +  ) 

and 

f  3  .  2  fQ2 

C2  =  $“U<X2  2^a2  +  ~  ^ 

Using  v2x  »  D  and  assuming  a  solution  form  for  x  as 

x  a  cosa2x( CjcosZaiy  +  C4)  (4 

then  C3  and  C4  become 

C3  =  -  Ci/(4ai2  +  a22) 

and 

C4  3  -  C2^aZ^ 
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Observe  that  <j>  equals  a  constant  and  equals  zero  along 
the  north/south  walls.  With  the  expressions  for  stream- 
function,  (4-32).  and  velocity  potential,  (4-36),  appropri  ate 
formulae  for  u  and  v  can  be  derived. 


u  *  U  -  sin  <*2x(Aai  sin2«iy  +  02(03  cos  2ojy  +  C4))  (4-37) 


and 


cos  °2X( °22’( l"COS  2a^y)  -  2C3a^sin  2o-jy) 


(4-38) 


and  <t>  is  given  by  equations  (4-31)  and  (4-32). 

2.  Source  Term 

The  addition  of  a  source  term  into  the  continuity 
equation  will  test  the  resolvability  of  the  GFEM  model  for  a 
source.  The  source  is  constructed  to  represent  a  small- 
scale  meteorological  phenomenon  to  provide  a  measure  of  the 
response  of  the  GFEM  model  to  forcing  near  the  smallest  grid 
scale.  The  source  initial  conditions  must  be  able  to 
describe  a  small-scale  feature  embedded  in  a  large  synoptic 
situation.  The  initial  conditions  are  derived  by  choosing  a 
desired  solution  and  then  back  substituting  to  derive  the 
form  of  the  source  required.  This  is  a  unique  approach  when 
adding  a  source  term.  First,  a  source,  S,  is  added  to  the 
continuity  equation  in  the  form  of 


3<fr 

at 


+ 


M 

ay 


a(—  +  ^-) 
?v3x  3y' 


3  V 


S  $ 


(4-39) 
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Following  the  scale  analysis  of  Chapter  3-2  In  Haltlner  and 
Williams  (1980), a  new  quasi -geostrophl c  potential  vortlclty 
equation  Is  formed  with  the  source  term  Included, 


u  —  +  w  —)  l  r  -  +  f  s 

%  3x  *  ayMC  T  a 


0 


Solving  for  S  leaves 

s-  <£-!?£ +  A>>  ‘4-40] 

x  y 

The  technique  for  testing  the  source  term  is  to  pick  a 
streamf uncti on  and  then  solve  for  the  source  term  which  will 
satisfy  that  streamfunction.  The  particular  streamf uncti on 
evaluated  will  be  a  steady  state  solution.  The  source  in 
this  case  is 

s  =  _L(i±  JL/ill  .  i£i)  -  14  +  1^4))  it  ai 

Vs*  sx  »2  3y2’  9X  9Ax2  3y2  (4‘4’ 

The  particular  form  of  the  streamfunction  is 

2  -  sin2^) 

*  3  -  U(y-ymid)  -  Asin  ^  e  R  L  (4-42 

where  R  is  a  constant.  The  use  of  the  exponential  allows  a 
variation  in  the  scale  whereby  a  small  scale  feature  can  be 
embedded  In  a  larger  flow  pattern.  Substitution  of  this 
particular  y  into  the  expression  for  the  source  yields 


•  V1"2  ¥ 
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E.  NONLINEAR  EFFECTS 


Isin 


(4-43) 


Cullen  (1974a)  discussed  the  requirement  to  simulate 
the  cascade  of  energy  to  unresolved  wavelengths.  To  avoid 
accumulation  of  energy  in  the  shortest  resolvable  wave¬ 
lengths,  he  developed  a  two-step  method  for  computing  the 
nonlinear  terms.  Implementation  of  the  two-step  method  can 
be  accomplished  efficiently.  In  the  formal  development  of 
the  model  equations,  there  are  several  nonlinear  terms.  The 
two-step  method  involves  projecting  both  variables  in  a 
nonlinear  term  into  Galerkin  space  and  then  performing  the 


multiplication.  This  procedure  insures  the  best  approxima¬ 
tion  to  that  nonlinear  term  in  a  weighted-mean  sense. 
Additionally,  the  computational  effort  for  the  two-step 
method  is  much  less  than  the  one-step  method  which  involves 
computation  of  all  the  possible  interaction  coefficients. 

F.  STABILITY  CRITERION 

The  stability  criterion  for  the  forecast  equations  is 


This  is  more  restrictive  than  a  normal  Courant-Fredri ch-Lewy 
( CFL)  condition.  It  is  made  up  of  the  normal  CFL  criterion 
based  on  advection  plus  another  effect  due  to  rotation 
(inertial  motion).  It  provides  an  upper  boundary  to  the 
maximum  allowable  time  step  when  using  a  semi -i  mpl  i  ci  t 
scheme.  In  this  research,  time  steps  in  excess  of  an  hour 
are  easily  used.  Care  must  be  taken  not  to  exceed  the 
stability  criterion  when  larger  mean  flows  are  experienced 
or  when  moving  to  higher  latitudes  (larger  f ) . 

A  more  detailed  discussion  of  the  derivation  of  the 
stability  criterion  can  be  found  in  Appendix  A. 
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EXPERIMENTS  AND  RESULTS 


Experiments  will  be  performed  Involving  the  features  in 
Chapter  II  to  test  the  hypothesis.  Four  separate  experi¬ 
ments  will  be  performed.  The  first  experiment  will 
establish  the  overall  performance  characteristics  of  three 
specified  models.  The  next  three  experiments  will  address 
each  specific  feature  supporting  the  hypothesis. 

Standards  of  comparison  must  be  established  for  each 
experiment.  The  basic  comparison  will  be  achieved  through 
harmonic  analysis  of  the  initial  and  forecast  fields.  The 
harmonic  analysis  is  a  double  Fourier  decomposition.  The 
harmonic  analysis  requires  that  the  data  be  on  a  uniform 
grid.  Because  of  the  triangular  and  variable  grids,  an 

a» 

interpolation  is  necessary.  A  fifth  degree  polynomial 

surface  fitting  routine  from  the  International  Mathematics 

and  Statistics  Library  (IMSL)  was  employed  to  interpolate 

data  on  the  non-uniform  grids  to  a  uniform  spacing.  First, 

a  harmonic  analysis  is  performed  along  columns  of  constant  x 

to  obtain  the  y  structure  of  the  initial  conditions.  Then  a 

harmonic  analysis  of  the  amplitudes  of  each  y  wave  number 

along  rows  gives  a  true  double  harmonic  analysis.  This 

double  transform  approach  is  necessary  due  to  the  initial 

2 

conditions.  The  use  of  sin  in  the  y  direction  is  the 
same  as  using  (1.0  -  cos  -jp).  The  boundary  conditions  are 

\ 


I 


satisfied  by  a  sine  wave  in  the  y  direction.  Therefore,  an 
infinite  sine  series  is  needed  to  represent  the  cosine 
functi on . 
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The  theoretical  movement  of  each  wave  is  obtained  from 
an  analysis  of  the  shallow-water  equations.  The  analysis  in 
Section  2-6  of  Haltiner  and  Williams  (1980)  predicts  the 
wave  speed  as 


C 


U  + 


(f/H)|| 
mZ  +  T/gH 


(5-1) 


where 

C  equals  wave  speed, 

H  equals  height,  and 
u  x  direction  wave  number. 

By  generalizing  the  streamf uncti on  to  a  two-dimensional  wave 
form 

4;  =  A  sin  cos  u(x-ct)  (5-2) 

and  after  further  manipulation,  the  wave  speed  becomes 

c  *  uci/Ci  +  f02/KMx2  +  uy2))) 

where 

$  equals  average  geopotential  height, 
ux  equals  x  direction  wave  number,  and 
My  equals  y  direction  wave  number. 
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A.  EXPERIMENT  1 

Establishing  wave  propagation  characteristics  cn  a 
uniform  grid  will  illuminate  the  properties  of  each  model. 
Three  models  will  be  compared: 

1.  A  triangular  subdivision  GFEM  model; 

2.  A  rectangular  subdivision  GFEM  model;  and 

3.  An  equivalent  finite  difference  model. 

The  GFEM  models  will  use  the  differentiated  form  of  the 
shallow-water  equations  with  a  semi -impl i ci t  scheme.  The 
finite  difference  model  uses  the  undifferentiated  form  on  a 
staggered  grid  with  a  centered  (leapfrog)  time  scheme.  The 
triangular  subdivision  uses  equilateral  triangles.  The 
relationship  between  base  and  height  for  an  equilateral 
triangle  is 

BASE  *  2  x  HEIGHT//! 

The  x  distance  therefore  has  the  same  relationship  to  the  y 
di stance 

ax  =  2  ay//3 

To  perform  comparisons  of  similar  models,  the  same  x,  y 
relationship  is  retained  for  the  rectangular  and  finite 
difference  models. 

The  basic  difference  between  the  rectangular  and 
triangular  models  will  be  in  the  approximating  polynomials. 
The  rectangular  polynomials  are  bilinear  while  the 
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triangular  polynomials  are  linear.  The  higher  order 
polynomials  should  provide  an  increase  in  accuracy.  During 
the  integration  process,  many  integrals  require  evaluation. 
Numerical  quadrature  could  be  utilized,  however,  a  more 
efficient  method  is  available  through  the  use  of  natural 
coordinates.  Quadrature  with  no  error  can  be  accomplished 
by  formula  with  this  method.  A  detailed  description  of  the 
natural  coordinate  method  is  delineated  in  Appendix  B  for 
both  triangles  and  rectangles. 

Diagrams  for  the  grids  for  the  triangular,  rectangular 
and  finite  difference  models  are  shown  in  Figures  2 ,  3  and  4 
respectively.  The  domain  is  4896.0  km  in  the  y  direction 
and  5653.0  km  in  the  x  direction.  Each  model  has  12  incre¬ 
ments  in  the  x  and  y  directions.  This  gives  156  degrees  of 
freedom.  All  tests  will  be  for  a  48  hour  time  integration, 
a  mean  flow  of  10  m/s  and  a  mean  depth  of  1000  meters.  A 
perturbation  of  1  m/s  will  be  added  to  the  geopotential 
field,  which  includes  the  mean  height  plus  the  north/south 
slope  required  for  the  10  m/s  mean  flow.  This  small  pertur¬ 
bation  will  focus  on  the  linear  aspects  of  the  model's 
capabilities.  Experiments  described  below  will  evaluate 
larger  perturbations,  and  hence  the  nonlinear  interactions. 
Models  are  evaluated  for  wave  numbers  1,  2,  3  and  4. 
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B.  RESULTS  1 

The  propagation  characteri sties  of  the  GFEM  models  and 
the  finite  difference  model  on  a  uniform  grid  will  be 
evaluated.  Graphical  displays  will  highlight  synoptic 
characteristics.  Harmonic  analysis  will  show  phase 
propagation  and  wave  amplitude  distributions. 

The  initial  and  48-h  forecasts  for  the  conditions 
stated  in  experiment  1  with  wave  number  one  on  a  triangular 
grid  are  shown  in  Figures  5  and  6.  The  southeast/northwest 
orientation  in  the  forecast  fields  is  due  to  the  sine  wave¬ 
form  in  the  y  direction.  This  waveform  with  model  boundary 
conditions  requires  an  infinite  number  of  y  modes. 
Fortunately,  the  amplitudes  of  the  higher  wave  numbers  are 
small.  Consequently,  there  are  only  a  few  important  modes. 
The  multiple  modes  cause  a  skewness  in  the  solutions  for  all 
three  models  since  the  y-modes  have  different  phase  speeds. 
Additionally,  an  analytic  expression  for  the  solution 
requires  an  infinite  sum  of  the  modes.  For  this  reason,  no 
true  solution  has  been  computed  to  compare  with  the  GFEM  and 
FDM  models. 

The  divergence  equation  is  the  most  sensitive  equation 
of  the  three  forecast  equations.  A  mean  depth  of  1000 
meters  was  chosen  because  the  divergence  is  large.  There¬ 
fore,  a  critical  evaluation  of  the  sensitivity  of  the  GFEM 
model  to  this  Important  meteorological  parameter  can  be 
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Initial  conditions  for  the  GFEM  model  with  a 
triangular  subdivision  and  wave  number  one. 
Contour  intervals  are  600  m2/s2  for  geopotential 
height,  .2  m/s  for  u  and  v,  .6x10"°  s"*  for 
vorticity  and  .2xl0"'s"^  for  divergence.  Nodal 
points  are  denoted  by  an  x. 


performed.  A  divergence  field  which  has  propagated 
uniformly  with  the  vorticity  field  is  shown  in  Figure  6. 
Harmonic  analysis  (not  shown)  indicates  divergence  ampli¬ 
tudes  which  initially  adjust  and  then  oscillate  slightly 
about  a  mean  value  of  .1270x10“^  s"*.  The  phases  of  all 
forecast  fields  are  uniform.  Synopti cal ly,  the  forecast 
fields  4>,  u,  and  v  indicate  a  smoothly  resolved  atmospheric 
wave . 

A  rectangular  subdivision  48-h  forecast  for  wave  number 
one  is  shown  in  Figure  7.  Comparisons  between  Figures  6  and 
7  show  little  if  any  difference  which  indicates  that  the 
triangular  and  rectangular  subdivisions  give  comparable 
results.  Harmonic  analysis  of  the  divergence  fields 
indicates  a  small  oscillation  about  a  mean  similar  to  that 
observed  for  the  triangular  subdivision.  The  48-h  forecast 
from  the  finite  difference  model  for  wave  number  one  is 
shown  in  Figure  8.  The  equivalent  finite  difference  model 
also  has  12  increments  in  each  direction.  Comparisons 
between  Figures  6,  7  and  8  show  that  the  finite  difference 
model  vorticity  field  lags  the  GFEM  models.  Additionally, 
the  divergence  field  has  lost  its  areal  extent  and  developed 
a  strong  gradient  of  divergence. 

In  each  of  the  wave  number  one  tests  there  were  12  grid 
points  representing  the  wave.  Good  propagation  should  be 
expected  from  all  three  of  the  models  for  a  12-increment 


wave.  As  the  wave  number  is  Increased  with  the  same  degrees 
of  freedom,  the  wave  resolution  will  be  decreased,  and 
clearer  comparisons  can  be  made. 

Figures  9  and  10  are  the  initial  and  48-h  GFEM  model 
forecasts  for  wave  number  2,  triangular  subdivision. 
Forecasts  for  the  rectangular  subdivision  and  the  FDM  model 
are  shown  In  Figures  11  and  12.  The  triangular  and 
rectangular  subdivision  forecasts  are  similar.  However,  the 
FDM  model  forecast  displays  additional  divergence  and 
vortlcity  centers  near  the  boundaries.  Figures  13  through 
16  are  a  series  similar  to  Figures  9  through  12  except  for 
wave  number  three.  Thece  waves  are  represented  by  4  grid 
points.  There  are  indications  in  Figure  14  that  energy  is 
appearing  in  wave  numbers  other  than  wave  number  3. 
Comparison  of  Figures  14  and  15  show  that  the  triangular 
forecast  contains  more  noise  than  the  rectangular  forecast. 
Figure  16  shows  a  divergence  field  from  the  FDM  model  which 
Is  very  different  from  Figures  14  and  15.  There  is  a 
north/south  asymmetry  in  the  u  and  divergence  fields.  The  u 
cells  In  the  northern  regions  appear  to  be  expanding  while 
the  southern  cells  are  decreasing. 

Table  lisa  composite  of  the  harmonic  analysis  for 
the  triangular,  rectangular  and  finite  difference  wave 
number  3  case  (v  component  amplitudes).  The  Initial  and 
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Fig.  14.  As  in  Fig.  13  but  a  48-h  forecast 
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48-h  v  component  amplitudes  for  wave  number  3  in  the  x 
direction  and  the  first  six  possible  y  modes  are  tabulated 
bel ow. 

TABLE  1 


Harmonic  analysis  (v  component  m/s)  -  wave  number  3 


Tri angl es 

Rectangl es 

Finite  D  i 

i  f ference 

Wave 

Initial 

48-H 

Initial 

48-H 

Initial 

48-H 

1 

2.12 

2.14 

2.20 

2.22 

2.20 

2.20 

2 

.01 

.04 

.01 

.02 

.01 

.09 

3 

.41 

.43 

.44 

.44 

.44 

.42 

4 

.00 

.00 

.00 

.00 

.00 

.00 

5 

.05 

.06 

.06 

.06 

.06 

.05 

6 

.00 

.00 

.00 

.00 

.00 

.00 

The  minor  differences  in  the  initial  state  amplitudes 
are  due  to  the  differences  between  the  triangular  and 
rectangular  approximation  for  the  initial  conditions. 
However,  the  48-h  amplitudes  show  that  the  rectangular 
version  has  slightly  better  preserved  the  initial  state 
amplitudes  than  either  the  triangular  or  finite  difference 
models.  The  finite  difference  model  has  transferred  more 
amplitude  to  the  other  modes,  especially  in  the  second  y 
mode. 

Figures  17  and  18  are  the  initial  and  48-h  GFEM  model 
triangular  subdivision  forecasts  for  wave  number  4.  Wave 
number  4  is  a  severe  test  for  this  model.  The  energy 
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g.  18.  As  in  Fig.  17  but  a  48-h  forecast. 
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transfer  to  other  modes  Is  readily  apparent.  Additionally, 
both  the  u  component  of  the  Mind  and  the  divergence  show 
considerable  noise  at  the  boundaries.  Figure  19  is  the  48-h 
forecast  for  the  rectangular  model.  No  energy  transfer  nor 
noise  at  the  boundaries  is  apparent  in  Figure  19.  Figure  20 
is  the  48-h  forecasts  for  the  finite  difference  model.  The 
asymmetry  in  the  wave  number  three  finite  difference  case 
is  magnified  for  wave  number  four.  The  divergence  field  is 
heavily  asymmetric.  The  u  cells  have  developed  a  dipole 
from  each  original  center.  Comparisons  of  Figures  18,  19 
and  20  show  the  superiority  of  the  rectangular  forecast. 

While  the  triangular  version  is  noisy,  it  still  has 
maintained  the  main  features  better  than  the  finite 
di f ference  model . 

Table  2  is  similar  to  Table  1  except  for  the  wave  number 
4  case. 

TABLE  2 

Harmonic  analysis  (v  component  m/s)  -  wave  number  4 

Triangles  Rectangles  Finite  Difference 


Wave 

Initial 

48-H 

Initial 

48-H 

Initial 

48-H 

1 

2.68 

2.60 

2.94 

2.96 

2.94 

2.96 

2 

.00 

.03 

.01 

.10 

.01 

.05 

3 

.51 

.59 

.58 

.56 

.58 

.60 

4 

.00 

.00 

.00 

.04 

.00 

.05 

5 

.06 

.08 

.08 

.08 

.08 

.08 

6 

.00 

.00 

.00 

.00 

.00 

.01 
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Fig.  20.  As  in  Fig.  18  but  a  FDM  model 


The  harmonic  analysis  shows  that  the  triangular  version 
transfers  more  energy  in  the  fifth  y  mode  than  does  either 
other  model.  However,  the  rectangular  model  transfers  more 
energy  in  the  second  y  mode  than  does  either  other  model. 
Table  3  compares  the  divergence  for  the  wave  number  4  case. 

TABLE  3 

Harmonic  analysis  (divergence  s-*)  -  wave  number  4 
(All  numbers  are  scaled  10  to  the  minus  10.) 

Triangles  Rectangles  Finite  Difference 


Wave 

Initial 

48-H 

Initial 

48-H 

Initial 

48-H 

1 

2682.0 

1741.0 

2929.0 

2179.0 

2801.0 

4462.0 

2 

.0 

140.8 

.0 

207.5 

447  .4 

5934.0 

3 

521.3 

391.1 

589.8 

440.7 

545.5 

2314.0 

4 

.0 

20.5 

.0 

81.2 

.0 

923.5 

5 

62.6 

60.5 

84.4 

56.8 

68.7 

1382.0 

6 

.0 

10.5 

.0 

8.4 

.0 

882.0 

The  finite  difference  model  has  greatly  amplified  the 
initial  divergence.  A  choice  between  the  triangular  or 
rectangular  grids  version  based  only  on  these  results  would 
be  purely  subjective.  However,  the  response  of  both  grids 
at  this  high  wave  number  in  terms  of  divergence  is  highly 
superior  to  the  FDM. 
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Figure  21  Is  the  phase  propagation  dla^am  for  the 
three  models  obtained  from  the  double  harmo.jc  analysis.  In 
all  cases,  the  percentage  propagation  Is  for  the  first  y 
mode  and  the  appropriate  x  mode.  The  phase  propagation 
diagram  Indicates  that  the  triangular  and  rectangular  models 
are  comparable  for  the  low  wave  numbers  and  diverge  slightly 
at  higher  wave  numbers.  In  all  cases  the  GFEM  models  have 
better  wave  propagation  than  the  FDM  model. 

The  analysis  performed  on  the  uniform  grids  has  shown 
that  each  GFEM  model  performs  in  a  highly  satisfactory 
fashion.  There  are  minor  differences  in  the  phase 
propagation,  while  the  rectangular  version  has  an  advantage 
in  the  control  of  energy  transfers  which  manifest  themselves 
in  the  divergence.  The  finite  difference  model  performs  as 
well  as  the  GFEM  models  for  the  long  waves.  However,  when 
forecasting  the  shorter  waves,  the  finite  difference  model 
does  not  perform  as  well  as  either  GFEM  model. 

C.  EXPERIMENT  2 

This  experiment  will  investigate  two  options  for  a 
variable  grid.  Atmospheric  wave  propagation  on  variable 
triangular  and  rectangular  grids  will  be  the  test  vehicle. 
Since  this  FDM  model  has  no  capability  on  a  variable  grid, 
it  will  not  be  evaluated  during  this  experiment.  The 
comparisons  will  be  between  a  rectangular  and  triangular 
GFEM  model . 
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Phase  propagation  diagram  for  the  triangular, 
rectangular  and  finite  difference  models  from 
Experiment  1. 


There  are  many  criteria  for  determination  of  a  suitable 
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variable  grid.  Both  rectangles  and  triangles  allow  a 
successful  implementation  of  a  variable  grid.  One 
criterion  is  the  accuracy  of  resolving  the  atmospheric  waves 
that  move  through  the  variable  grid.  A  second  criterion  is 
the  increase  in  resolution.  Another  consideration  is  the 
ease  with  which  the  grid  can  be  refined. 

Cullen  (1979)  stressed  the  importance  of  a  smooth 
variation  in  element  size  when  increasing  the  resolution 
from  coarse  to  fine.  Older  (1981)  developed  a  technique  for 
transforming  a  uniform  triangular  grid  into  a  variable  grid 
when  working  with  a  GFEM  model.  He  demonstrated  that  a 
smoothly  varying  grid  allowed  atmospheric  waves  to  propagate 
through  the  transition  zone  with  much  less  noise  generation 
than  an  abruptly  varying  grid.  Some  finite  difference 
models  with  nested  grids,  for  example  Harrison  (1973),  have 
abruptly  varying  resolution.  A  general i zati on  of  Older's 
techniques  allows  a  similar  transformation  for  a  rectangular 
gri  d . 

The  test  cases  for  experiment  2  will  be  for  a  simple 
atmospheric  wave.  A  mean  depth  of  5  km,  a  mean  flow  of 
10  m/s,  wave  number  one  and  a  perturbation  of  1.0  m/s  are 
empl  oyed . 
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D.  RESULTS  2 

The  ability  to  refine  a  grid  should  be  easily  achieved. 
The  method  of  transforming  from  a  uniform  grid  to  smoothly 
varying  should  allow  one  to  choose  not  only  what  degree  of 
refinement  but  also  where  the  refinement  occurs.  Older 
(1981)  developed  a  transformation  method  based  on 

X  *  x  +  A  cos  Kx 


where 

X  3  the  transformed  grid 
x  3  the  original  grid 
A  3  a  constant 
k  3  2w/L. 

Woodward  (1981)  modified  the  transf ormati on  to  include  a 


for  the  longitudinal  stretching.  A  trignometric  identity 
is  employed  yielding 

X  =  x  +  A(l  -  cos  kx) 

This  research  added  the  capability  to  control  the  location 

of  the  refinement  through 

X  3  x  +  A(1  -  cos(kx  +  5)) 
aX 

The  map  factor  is  defined  by 

|4  =  1  +  kA  sin(kx  +  6) 
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The  maximum  and  minimum  values  are 


3X 

3X. 


■  1  +  kA,  f£  -  1  -  kA 


"max  “min 

The  ratio  R  of  maximum  map  factor  to  minimum  map  factor  is 


R 


1+kA 

’ra 


Sol vi ng  for  A  yields 

.  R-l 
A  =  k(R+l ) 

For  purposes  of  this  experiment  R  was  chosen  to  vary  from 
1.0  to  4.0.  A  value  of  R  *  4  implies  that  the  minimum  X  is 
one  fourth  of  the  maximum  X.  The  Y  transformation  is 
performed  in  a  similar  fashion  except  that  Y  =  y  +  B  sin  Ly 
is  employed  where 

Y  »  the  transformed  grid, 
y  =  the  original  grid, 

B  *  a  constant,  and 
l  *  2*/W 

Placement  of  the  high  resolution  is  accomplished  by 
determining  the  value  of  6  required  to  place  the  minimum  map 
factor  as  desired.  For  instance,  if  the  refined  grid  is 
desired  in  the  middle  of  the  channel  (L/2)  then  5  would  be 
equal  to  x/2. 
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Figure  22  is  the  48-h  forecast  with  the  triangular 
subdivision  for  a  uniform  grid.  This  forecast  will  act  as 
the  control  for  the  triangular  cases.  Figure  23  is  the 
triangular  subdivision,  48-h  forecast  with  R  *  2.0.  Figure 
24  is  the  triangular  subdivision,  48-h  forecast  with 
R  *  3.0.  Figure  25  is  the  triangular  subdivision,  48-h 
forecast  with  R  *  4.0.  All  of  the  48-h  u  and  v  fields 
appear  identical.  This  Indicates  that  there  is  no  degra¬ 
dation  in  forecast  skill  when  using  the  same  number  of 
degrees  of  freedom  and  locating  many  of  those  unknowns  into 
a  refined  area.  There  is  no  noise  apparent  in  the 
transition  regime.  However,  there  are  definite  differences 
in  the  divergence  field,  although  the  magnitude  of  the 
maximum/  minimum  value  of  divergence  is  relatively  small 
( 1 .0 xl 0  8  s'1)- 

Table  4  shows  the  harmonic  analysis  amplitudes  cf  the 
geopotential  fields  for  the  different  triangular  cases.  The 
relationship  between  initial  and  final  amplitudes  remains 
the  same  regardless  of  the  degree  of  variability.  The  phase 
propagation  of  the  v  field  in  the  control  case  is  97.7%, 
whereas  it  is  96.7%,  95.4%  and  94.1%  for  R  values  of  2,  3 
and  4.  The  variation  from  the  control  of  the  phase  propaga¬ 
tion  with  grid  refinement  is  under  four  percent. 
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Fig.  24.  As  in  Fig.  22  for  R  *  3.0, 
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TABLE  4 


Harmonic  Analysis  (Geopotential  m^/s^) 
Tri angles 


R  * 

1.0 

R  * 

2.0 

R  * 

3.0 

R  = 

4.0 

Wave 

Initial 

48-H 

Initial 

48-H 

Initial 

48-H 

Initial 

48-H 

1 

67.98 

67.16 

67.91 

64.65 

67.91 

64.95 

67.91 

63.70 

2 

.00 

.55 

.00 

.69 

.00 

.99 

.00 

.57 

3 

13.55 

12.90 

13.45 

12.55 

13.41 

13.24 

13.29 

11.16 

4 

.00 

.21 

.00 

.37 

.00 

.54 

.00 

.28 

5 

1.90 

1.70 

1.67 

1.30 

1.33 

1.62 

1.32 

1.38 

6 

.00 

.16 

.00 

.23 

.00 

.36 

.00 

.21 

Figure  26  is  the  control  case  for  the  rectangular 
version,  48-h  forecast.  Figure  27  is  the  rectangular 
subdivision,  48-h  forecast  with  R  *  2.0.  Figure  28  is  the 
rectangular  subdivision,  48-h  forecast  with  R  *  3.0. 

Figure  29  is  the  rectangular  subdivision,  48-h  forecast 

with  R  =  4.0.  Close  inspection  reveals  no  noise  in  either 
the  t,  u  or  v  fields.  No  degradation  in  forecast  skill  has 
been  observed  due  to  increased  resolution  with  a  rectangular 
subdivision.  All  f,  u  and  v  fields  are  very  similar  in 
structure.  The  main  difference  in  the  forecasts  lies  in  the 
divergence  field  where  the  magnitudes  are  small  (  10"®s”M 
because  of  the  specified  mean  depth. 
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Fig.  27.  As  in  Fig.  26  for  R  »  2.0. 
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Table  5  is  similar  to  Table  4  except  for  the 
rectangular  cases.  Phase  propagation  for  the  control  case 
v  field  is  98.0%,  whereas  it  is  98.3%,  98.1%,  and  97.9%  for 
R  values  of  2,  3,  and  4.  The  variation  from  the  control 
case  of  the  phase  propagation  with  grid  refinement  is  less 
than  half  a  percent. 

TABLE  5 

Harmonic  Analysis  (Geopotential  m^/s^) 

Rectangl es 


R  = 

1.0 

R  = 

2.0 

R  = 

3.0 

R  * 

4.0 

Wave 

Ini ti al 

48-H 

Ini  ti al 

48-H 

Ini ti al 

48-H 

Ini ti al 

48-H 

1 

68.03 

67.48 

67.98 

65.26 

68.22 

66.45 

68.14 

66.34 

2 

.01 

1.00 

.03 

.36 

.00 

.85 

.00 

1.22 

3 

13.59 

12.95 

13.39 

13.26 

13.49 

13.81 

13.36 

12.39 

• 

4 

.01 

.46 

.05 

.22 

.01 

.48 

.00 

.43 

5 

1.92 

1.81 

1.57 

1.64 

1.54 

1.59 

1 .46 

1.19 

6 

.02 

.23 

.07 

.04 

.01 

.30 

.00 

.22 

The  tests  thus  far  indicate  a  successful  grid  refine¬ 
ment  capability  for  both  triangles  and  rectangles.  There 
are  neither  drastic  impairments  nor  improvements  in  either 
phase  propagation  or  amplitude  changes  for  either  technique. 
Neither  version  stands  out  above  the  other  although  the 
rectangular  version  has  less  phase  propagation  variation  and 
better  amplitude  conservation.  Both  versions  allow  compara¬ 
ble  grid  refinement  through  the  use  of  the  procedure 
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previously  described.  Both  allow  an  easy  and  straightfor¬ 
ward  implementation  of  the  procedure. 

However,  there  is  a  more  fundamental  difference  between 
the  versions.  Operators,  such  as  a  simple  derivative,  have 
weighting  coefficients  associated  with  finite  element  appli¬ 
cations  as  well  as  with  finite  difference  schemes.  For 
instance,  a  simple-centered,  one-dimensional  finite  differ¬ 
ence  derivative  has  weights  of  1/2  and  -1/2.  The  finite 
element  weights  are  identical  to  the  finite  difference 
weights  in  this  one-dimensional  case.  The  weights  can 
similarly  be  computed  for  two-dimensional  operators  such  as 
a  laplacian.  The  quadrature  formulas  used  to  determine  the 
weights  are  explained  in  more  depth  ii  Appendix  B.  The 
geometry  of  the  triangles  affects  the  weights.  For 
instance,  the  two-dimensional  Laplacian  applied  to  a  finite 
element  application  on  an  equilateral  triangular  subdivision 
gives  weights  of  1/6  to  all  surrounding  points,  and  -*  +o 
the  center  point,  as  indicated  in  Figure  30.  However,  the 
weights  change  if  a  triangle  with  a  height  equal  to  one  half 
the  base  is  employed  as  shown  in  Figure  31.  The  fact  that 
some  weights  are  zero  is  a  strong  cause  for  alarm.  If  the 
triangles  are  flattened  further,  the  farthest  end  point 
weights  become  negative.  This  is  strictly  a  function  of  the 
geometry  of  the  triangle  and  is  a  clear  sign  that  use  of 
the  triangles  warrants  extreme  caution.  In  fact,  Williams 


30.  Weights  associated  with  a  Laplacian  operator 

for  a  triangular  subdivision  using  equilateral 
tri angl es . 


1/4  1/4 


31.  Weights  associated  with  a  Laplacian  operator  for  a 
triangular  subdivision  using  triangles  with  a  base 
equal  to  twice  the  height. 


(unpublished  notes)  has  determined  that  the  geometry  of  the 
triangles  dictates  that  unless  a  triangle  with  at  least  two 
equal  sides  (isosceles)  is  employed  the  boundary  conditions 
will  not  be  exactly  satisfied. 

E.  EXPERIMENT  3 

The  capability  of  the  GFEM  model  to  resolve  atmospheric 
phenomena  near  the  scale  of  the  smallest  grid  length  will  be 
demonstrated.  This  experiment  is  in  response  to  Task  II  of 
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the  hypothesis. 

The  ultimate  motivation  for  any  "improvement"  in  a 
model  is  a  better  forecast.  During  the  past  decade, 
improvements  have  come,  to  a  certain  degree,  by  increased 
resolution.  While  it  is  intuitively  obvious  tnat  increased 
resolution  will  improve  a  forecast,  demonstrations  have  not 
been  made  to  show  how  the  model  resolves  features  near  the 
smallest  grid  length.  Staniforth  and  Mitchell  (1978) 
increased  the  resolution  with  a  variable  grid  but  resolved  a 
synoptic  scale  feature  in  a  37x37  uniform  fine  mesh  area. 

In  this  experiment,  a  source  term  simulating  a  mass 
source  is  added  to  the  continuity  equation.  An  expression 
for  the  source  term  in  terms  of  a  wave  form  is  derived  by 
theoretical  means.  The  wave  form  allows  a  very  small  scale 
feature  to  be  simulated.  A  detailed  description  of  this 
procedure  can  be  found  in  Chapter  IV-D-2.  By  assuming  the 
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wave  is  in  equilibrium  and  steady  state,  the  solution  and 
source  are  known  for  all  time  when  the  Rossby  number  is 
smal  1  . 

The  development  of  the  source  term  expression  followed 
traditional  quasi -geostrophi c  theory.  One  requirement  for 
quasi -geostrophi c  theory  is  that  the  Rossby  number  be  small. 
In  the  linear  case,  the  Rossby  number  will  be  approximately 
0.1.  However,  in  the  nonlinear  case  the  Rossby  number  is 
approximately  0.4  and  the  quasi -geostrophi c  theory 
assumption  is  violated. 

The  first  test  will  be  on  a  uniform  grid  for  the 
rectangular  GFEM  and  finite  difference  models  and  the  second 
test  wi'l  be  on  a  variable  rectangular  grid  GFEM  model. 

High  resolution  versions  of  the  uniform  test  will  be  run  and 
act  as  the  control.  A  mean  depth  of  5  km  was  selected. 
Perturbations  of  2.5  m/s  and  25  m/s  upon  a  mean  flow  of 
10  m/s  will  be  evaluated  for  the  uniform  grid.  This  will 
allow  both  the  linear  and  nonlinear  aspects  of  the  different 
models  to  be  observed.  Here  linearity  is  implied  by  the 
smallness  of  the  perturbation  amplitude.  The  source  term 
solution  is  nonlinear  anu  nonlinear  interactions  do  occur 
when  the  amplitude  is  small.  Each  test  will  be  integrated 
for  96  h  of  forecast  time.  A  perturbation  of  25  m/s  upon  a 
mean  flow  of  10  m/s  will  also  be  evaluated  for  the  variable 
grid. 
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The  main  thrust  of  this  experiment  is  to  determine  the 
grid-scale  sensitivity  for  the  GFEM  model  and  the  equivalent 
finite  difference  model.  Sensitivity  will  be  measured  by 
insertion  of  a  source  term  into  the  continuity  equation. 
However,  the  theoretical  development  presented  in  Chapter 
I V - D - 2  is  for  the  rotational  component  of  the  wind.  Because 
the  initial  conditions  for  u  and  v  are  non-di vergent ,  the 
divergence  in  the  rnodel  must  increase  from  a  zero  initial 
state.  This  increase  will  require  an  adjustment  process. 
This  serendipity  effect  will  be  exploited  during  the 
compari sons . 

The  linear  cases  examined  represent  a  small-scale  short 
wave  embedded  in  a  long  wave  pattern.  Using  a  small  per¬ 
turbation  for  the  first  tests  will  establish  a  level  of 
confidence  that  each  model  responds  to  the  source  term 
creditably.  The  nonlinear  cases  represent  an  active  vortex 
in  a  mean  flow.  The  maximum  u  component  is  35  m/s  which  is 
hurri cane/ typhoon  strength.  Physically,  the  source  term  in 
the  divergence  field  opposes  the  advection  of  the  vorticity 
field,  and  "anchors"  the  steady-state  vorticity  solution. 

In  either  the  linear  or  nonlinear  tests,  the  source 
term  is  nonlinear.  By  constraining  the  perturbation  to  be 
small,  the  source  term  plays  a  linear  role.  However,  when 
the  perturbation  is  large,  the  source  term  allows  full 
nonlinear  interactions. 
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F.  RESULTS  3 

1 .  Linear  Case  -  Unifora  Grids 

Figures  32  and  33  are  the  initial  and  96-h  fore¬ 
casts  for  the  rectangular  GFEM  model.  There  are  156  degrees 
of  freedom  in  these  forecasts.  The  initial  v  field  has  a 
maximum  value  of  +2.27  m/s  across  a  separation  of  4  incre¬ 
ments.  The  initial  u  field  maximum  and  minimum  are  12.5  m/s 
and  7.5  m/s.  The  source  term  is  shown  in  the  initial  fields 
and  remains  constant  throughout  the  integration.  The  diver¬ 
gence  value  increased  to  a  steady  state  solution  during  the 
first  12  hours  with  an  oscillation  that  died  out  after  hour 
24.  Figure  34  is  the  equivalent  finite  difference  96-h 
forecast.  Figure  34  demonstrates  that  the  finite  difference 
model  performs  well  for  these  initial  conditions.  Figure  35 
is  the  high-resolution  96-h  forecast  for  the  rectangular 
GFEM  model.  Close  inspection  of  Figures  33  and  35  shows 
little  difference.  The  low  resolution  forecast  with  156 
degrees  of  freedom  has  converged  to  the  high  resolution 
forecast  with  600  degrees  of  freedom.  Figure  36  is  a  graph 
of  v  component  amplitudes  at  hour  96  as  a  function  of  x  and 
y  wave  numbers  for  the  low  and  high  resolution  and  finite 
difference  models.  The  high  resolution  (S253115)  and  low 
resolution  ( SI 531 15)  forecasts  are  In  excellent  agreement. 
However,  the  FDM  (F153115)  forecast  departs  from  the  control 
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Fig.  32.  Initial  conditions  for  the  GFEM  model  with  a 

rectangular  subdivision  and  source  term  added  to 
the  continuity  equation.  Resolution  Is  low  and 
uniform  for  a  2.5  mis  oerturbatl on.  Contour 
Intervals  are  600  m2/ sz  for  geopotentlal  height, 
.5  a/s.for  u  and  v,  .2x10“^  s'*  f°r  vortlclty  and 
.2x10“®  s“*  for  source. 
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Fig.  33.  As  In  Fig.  32  for  a  96-h  forecast. 
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j  Fig.  36.  V  component  amplitude  as  a  function  of  x  and  y 

!  wave  number  at  hour  96  for  the  high  (S253115)  and 

low  ( S1531 IS)  resolution  6FEM  models  and  the 
finite  difference  (F153115)  model  and 
1  perturbation  «  2.5  m/s.  Contour  Interval  Is 
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even  with  this  relatively  linear  case.  For  unlforalty  of 
comparison,  similar  ranges  of  x,  y  wave  numbers  will  be 
shown  In  Figures  36  and  41. 

2.  Nonlinear  Case  -  Variable  6r1ds 

Figures  37  and  38  are  the  Initial  and  96-h 
forecasts  for  the  rectangular  6FEN  model  with  low 
resolution.  The  Initial  v  field  has  a  maximum  of  ±22.7  a/s 
across  a  separation  of  4  Increments,  but  most  of  the  change 
occurs  over  2  Increments.  The  Initial  u  field  maximum/ 
minimum  Is  35  m/s  and  <15  m/s.  The  source  term  as  shown 
remains  constant  during  the  Integrations.  The  adjustment 
process  Is  evident  when  viewing  96-h  forecasts.  The  final  v 
field  has  a  maximum  of  44.3  m/s,  and  a  minimum  of  -38.0  m/s. 
The  final  u  field  has  a  maximum  of  40.2  m/s  and  a  minimum  of 
-29.8  m/s.  Figure  39  Is  the  equivalent  finite  difference 
model  forecast  of  Figure  38  and  the  marked  differences  In 
the  forecast  are  readily  apparent.  The  maximum  final  v 
component  Is  41.2  m/s  and  the  minimum  is  -47.3  m/s.  The 
maximum  final  u  component  Is  50.6  m/s  and  the  minimum  value 
is  -38.6  m/s.  Figure  40  Is  the  96-h  forecast  for  the 
rectangular  high-resolution  6FEM  model.  Comparisons  of 
Figures  38,  39  and  40  show  a  convergence  of  the  low- 
resolution  GFEM  version  towards  the  high-resolution  solution 
with  a  much  poorer  showing  for  the  finite  difference  model. 
Figure  41  Is  a  graph  of  v  component  amplitudes  at  hour  96  as 
a  function  of  x  and  y  wave  numbers  for  the  low  and  high 
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Fig.  41.  V  component  amplitude  as  a  function  of  x  and  y 

wave  number  at  hour  96  for  the  high  (S258115)  and 
low  (R1S8115)  resolution  GFEM  models  and  the 
finite  difference  (F158115)  model.  Perturba¬ 
tion  »  25.0  m/s.  Contour  Interval  Is  0.5  m/s. 
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resolution  and  finite  difference  model.  As  In  the  linear 
case,  there  is  strong  agreement  between  the  low  and  high 
resolution  forecasts.  The  departure  of  the  FOM  from  the 
control  Is  readily  apparent  in  this  highly  active  case.  The 
inability  of  the  FDM  to  properly  handle  the  nonlinear 
interactions  manifest  itself  as  large  spurious  amplitudes  at 
high  x  and  y  wave  numbers. 

The  relatively  poor  showing  of  the  finite  difference 
model  for  the  forced  cases  shown  so  far  can  be  Improved,  to 
a  certain  degree,  by  increasing  the  resolution.  However, 
increased  resolution  requires  more  computational  effort. 
Figure  42  is  the  96-h  forecast  with  a  perturbation  of  25.0 
m/s  for  the  FDM  model.  The  initial  conditions  are  identical 
to  the  nonlinear  cases  shown  in  Figure  37.  The  difference 
in  the  forecasts  is  the  resolution  where  there  are  now  576 
degrees  of  freedom  (24X24).  The  96-h  forecast  is  certainly 
better  than  the  12X12  forecast.  However,  there  is  some  high 
frequency  noise  in  the  forecast.  Figure  43  is  identical  to 
Figure  42  except  that  there  are  1296  degrees  of  freedom 
(36X36).  The  high  frequency  noise  is  readily  apparent  and 
is  a  manifestation  of  nonlinear  aliasing.  For  this  forced 
case  the  finite  difference  model  is  unable  to  achieve  the 


same  forecast  as  the  low  resolution  GFEM  model.  This  result 
will  be  amplified  when  a  model  with  t^e  same  number  of 
degrees  of  freedom  as  the  low  resolution  model  Is  employed 
but  with  a  variable  grid. 
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This  comparison  of  the  finite  element  and  finite 
difference  methods  for  the  source  term  dramatically  high- 
lights  the  potential  Improvement  when  utilizing  a  Galerkin 
approach.  Both  graphical  displays  and  harmonic  analysis 
leave  little  doubt  as  to  the  superiority  of  the  GFEM  model 
over  the  equivalent  finite  difference  model. 

Adding  an  extra  margin  to  the  already  established 
superiority  can  be  accomplished  by  using  the  low-resolution 
model  with  a  variable  grid.  Figure  44  Is  similar  to  the 
tests  In  Figure  38  except  that  a  variable  resolution  grid 
(R  *  2.0)  was  utilized.  A  96-h  forecast  for  a  variable 
resolution  grid  (R  *  3.0)  Is  shown  In  Figure  45.  Comparison 
of  Figures  38,  40,  44,  and  45  show  the  convergence  of  the 
low-resolution  uniform  solution  towards  the  high  resolution 
solution  as  the  degree  of  variability  is  Increased.  Figure 
46  1  s  a  graph  of  geopotential  amplitudes  at  hour  96  versus  y 
mode  wave  number  (x  wave  number  one)  for  the  uniform, 
variable,  and  high  resolution  models.  The  R  »  2.0  case 
better  represents  the  control  than  the  uniform  low 
resolution  forecast  especially  at  low  wave  numbers. 

However,  at  wave  number  three  and  above  all  models  have 
similar  differences  from  the  control. 

The  harmonic  analysis  demonstrates  an  Improvement  with 
Increasing  variability  but  Is  not  conclusive.  Another 
method  to  evaluate  the  Improvement  with  increased  varia¬ 
bility  Is  to  look  at  the  difference  charts  between  the 
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Fig.  46.  Geopotentlal  amplitude  versus  x  wave  number  1  and 
y  waves  1-6  at  hour  96  for  the  uniform,  varlabl 
and  high  resolution  GFEM  models.  Perturbation 
25.0  m/s. 
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control  and  each  level  of  variability.  The  difference  chart 
in  Figure  47  Is  the  control  96-h  forecast  (high  resolution) 
minus  the  low-resolution  96-h  forecast.  Figure  47  shows  a 
maximum  difference  In  the  center  of  the  vortex.  Figure  48 
Is  the  comparison  of  the  high  resolution  96-h  forecast  and 
the  R  «  2.0  variable  grid  96-h  forecast.  The  maximum 
difference  has  decreased  from  the  uniform  case.  Figure  49 
Is  the  comparison  of  the  high  resolution  96-h  forecast  and 
the  R  »  3.0  variable  grid  96-h  forecast.  The  difference 
centers  at  the  right-hand  edge  In  Figures  47,  48  and  49  are 
a  result  of  post  processing  Interpolation  errors.  Although 
the  maximum  difference  In  the  center  of  the  vortex  Is 
greater  In  Figure  49  than  Figure  48,  the  latitudinal 
difference  gradient  Is  less.  Root  mean  square  differences 
from  the  high  resolution  control  are  652.2  m2/s2,  608.8 
m2/s2  and  469.7  m2/s2  for  the  low  resolution,  R  ■  2.0  and 
R  ■  3.0  tests  respectively.  The  forecasts  have  definitely 
been  improved  at  each  level  of  Increased  resolution. 

6.  EXPERIMENT  4 

This  experiment  will  demonstrate  a  moving  grid 
capability.  Harrison  (1973)  demonstrated  the  use  of  a 
moving  grid  In  his  finite  difference  model.  The  Inherent 
problems  with  moving  finite  difference  models  are  abrupt 
changes  In  resolution  causing  spurious  noise.  The  noise, 
unless  damped.  Is  magnified  If  left  unattended  every  time  a 
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Fig.  47.  Forecasts  and  difference  charts  at  hour  96  for  the 
uniform  high  (S258115)  and  low  (S158115)  resolu¬ 
tion  SFEM  models.  Perturbation  *  25.0  m/s. 

Contour  Interval  J  s  ,600  m*/s*  for  the  geopotentlal 
fields  and  200  mz/s*  for  the  difference  field. 
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Fig.  48.  As  In  Fig.  47  for  the  high  resolution  (S258115) 
and  the  variable  (S178115)  (R  •  2.0)  6FEM  node! 
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Fig.  49.  As  In  Fig.  47  for  the  high  resolution  (S258115) 
and  the  variable  (S198115)  (R  -  3.0)  SFEM  models 


grid  «oves.  In  a  GFEM  model  there  Is  a  smooth  transition 
Into  the  fine  grid  area  with  no  attendant  noise.  Therefore, 
moving  the  entire  grid  with  the  fine  grid  Is  simply  a  matter 
of  accurate  Interpolation. 

Three  tests  will  be  performed.  First,  uniform  low  and 
high  resolution  models  will  be  Integrated.  These  models 
will  not  use  a  moving  grid.  Then  a  model  with  the  same 
degrees  of  freedom  as  the  low-resolution  model,  but 
employing  a  moving,  variable  grid,  will  be  run.  The  high- 
resolution  model  will  act  as  the  control.  Identical  Initial 
conditions  will  be  used  for  all  three  tests.  These  condi¬ 
tions  are  similar  to  those  developed  for  experiment  3, 
except  that  the  source  term  will  be  zero  throughout  the 
course  of  the  Integration.  The  Initial  conditions  Include  a 
sharp  trough  with  a  relative  vortlclty  maximum  which  Is 
advected  by  the  mean  flow  In  the  absence  of  any  forcing  (the 
source) . 

Of  particular  concern  is  the  generation  of  noise 
Incurred  during  grid  movement.  Therefore,  a  forecast  of 
96-h  Is  made  during  which  this  grid  moved  eleven  times.  As 
previously  stated,  the  model  utilizes  absolutely  no  diffu¬ 
sion,  filters  or  friction.  If  spurious  modes  are  excited 
during  the  movement  of  the  grid,  they  will  appear  in  the 
harmonic  analysis. 
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The  procedure  of  moving  the  finite  element  grid  is 
achieved  by  displacing  all  the  nodes  one  unit  distance  in 
the  direction  desired.  This  retains  the  relative 
distribution  of  the  nodes  and  merely  moves  the  entire  grid 
as  an  entity.  Ramifications  of  this  approach  are: 

1.  The  areas  of  each  element  remain  the  same  and  the 
numerical  quadrature  coefficients  for  all  the 
associated  integrals  do  not  have  to  be  recomputed. 

2.  The  unit  distance  moves  should  be  the  smallest 
grid  increment;  and 

3.  Oirect  solvers  that  utilize  a  preprocessing 
procedure  based  on  the  grid  distribution  will  not 
require  the  preprocessi ng  to  be  recomputed  after 
each  move. 

The  most  important  aspect  of  moving  the  grid  is  to  accurate¬ 
ly  interpolate  the  new  field  values  for  the  new  grid  point 
position.  In  this  dissertation  a  standard  International 
Mathematics  and  Statistical  Library's  (IMSL)  bi-cubic  spline 
interpolating  routine  was  employed.  The  movement  will  occur 
only  in  the  x  direction  and,  therefore,  a  one-dimensional 
interpolator  suffices.  The  interpolator  accommodates  the 
periodic  boundary  condition.  However,  Interpolators  for 
non-periodic  boundary  conditions  are  available.  The  accu¬ 
racy  of  the  interpolator  is  exact  at  grid  points  for  a 
uniform  distribution.  Initial  experimentation  verified  that 
a  field  could  be  moved  with  no  loss  in  accuracy  for  a 
uniform  grid.  Experiment  4  will  test  its  capability  to  move 


H.  RESULTS  4 

Figures  50  and  51  are  the  initial  and  96-h  forecasts 
for  the  uniform  1 ow-resol uti on  model.  The  perturbation  is 
25.0  m/s,  mean  flow  10  m/s  and  mean  depth  5  km.  The  initial 
divergence  is  a  zero  field  but  spins  up  through  geostrophic 
adjustment.  The  northeast/ southwest  trough  orientation  is  a 
manifestation  of  the  nonlinear  interactions  occurring  for 
this  highly  active  perturbation.  Figure  52  is  the  96-h 
forecast  fo?  the  uniform  high-resolution  model.  The 
northeast/ southwest  trough  orientation  is  evident  in  this 
forecast  as  well.  Figure  53  is  the  96-h  forecast  for  the 
low-resolution,  moving  variable  grid  (R  »  2.0)  model. 

Time  steps  were  4320s,  3085.7s  and  3600s  for  the  low, 
high  and  variable  models,  respectively.  The  variable  model 
took  95  full  time  steps  and  two  hal  .*  steps  during  the  96-h 
forecast.  The  half  steps  are  required  only  to  start  the 
model  to  obtain  the  N  and  N+l  time  levels  in  the  leapfrog 
time  stepping.  The  grid  moved  11  times  during  the  forecast 
at  time  steps  8,  16,  24,  32,  40,  48,  56,  64,  72,  80  and  88. 

The  x  axis  origin  in  Figure  53  is  3530  km  indicating 
eastward  grid  movement.  The  vortex  is  still  in  the  fine 
mesh  area.  Additionally,  there  is  no  noise  in  the  forecast 
fields.  To  reiterate  a  previous  statement,  there  are  no 
filters,  diffusion  or  friction  terms  employed  during  this 
Integration.  Comparisons  among  Figures  51,  52  and  53  show 
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that  the  variable  grid  forecast  more  closely  approaches  the 
control  than  does  the  uniform  1 ow-resol utl on  model.  This  Is 
to  be  expected  since  the  fine  mesh  area  has  been  moved  to 
coincide  Ml th  the  active  region  of  the  forecast  domain. 

Table  6  shows  the  geopotentlal  harmonic  analysis  for 
the  three  models  after  a  96-h  forecast. 


TABLE  6 

Harmonic  Analysis  (Geopotentlal)  for  Experiment  4 

Units  are  m2/s2 


Wave 

Low  Resolution 

High  Resolution 

Moving  Grid 

1 

876.1 

904.7 

850.0 

2 

167.4 

152.9 

146.6 

3 

122.5 

176.9 

159.6 

4 

54.7 

39.0 

39.0 

5 

20.3 

10.2 

9.6 

6 

5.1 

7.6 

7.3 

Comparisons  of  the  higher  wave  number  amplitudes  show 
almost  identical  amplitudes  in  y  modes  four,  five  and  six 
for  the  moving  case  and  high-resol ution  uniform  case.  The 
amplitude  in  y  mode  six  indicates  that  the  movement  of  the 
grid  has  not  generated  unwanted  noise.  The  overall  improve¬ 
ment  for  the  moving  case  is  a  result  of  the  finer  grid's 
ability  to  better  resolve  the  small-scale  Interactions 


VI.  SUMMARY  AND  RECOMMENDATIONS 

Numerical  weather  prediction  requires  that  the 
equations  governing  atmospheric  flow  be  solved  numerically. 
These  flows  cover  a  large  range  of  scales.  Classical 
modeling  schemes  have  successfully  forecast  the  longer 
scales  and  Improvements  are  being  sought  for  the  smaller 
scales. 

In  this  dl ssertatlon,  the  GFEM  has  been  employed  for 
the  shal  low-water  equations  In  a  channel  domain.  The 
differentiated  form  of  the  equations  on  an  unstaggered 
finite  element  grid  with  a  semi -1 mpl 1 ci t  time  scheme  has 
been  utilized.  An  equivalent  finite  difference  model 
representing  the  conventional  approach  has  also  been 
utilized.  Analytic  Initial  conditions  representing  simple 
atmospheric  waves  and  a  more  complex  source  term  have  been 
used  to  test  the  different  models. 

The  hypothesis  of  this  dissertation  Is  that  the  GFEM  Is 
a  viable  option  for  numerical  weather  prediction.  The 
theoretical  foundation  of  the  GFEM  leads  to  a  minimization 
of  the  error  between  the  actual  equations  and  their  approxi¬ 
mation.  This  minimization  implies  a  basis  for  expected 
Improvement  In  forecast  capabilities.  Three  features  were 
addressed  to  support  that  hypothesis.  First,  alternatives 
for  variable  grids  were  Investigated.  Second,  forcing  near 


the  smallest  grid  scale  determined  If  the  GFEM  model 
resolved  the  smaller  scales  better  than  a  finite  difference 
counterpart.  Finally,  a  proof  of  concept  for  the  ability  to 
move  the  grid  was  demonstrated.  Each  task  adds  to  the 
usability  of  the  GFEM  and  Illuminates  the  well-rounded 
potential  that  it  contains. 

First,  the  ease  with  which  the  method  allows  variable 
grids  was  demonstrated.  In  all  cases,  no  numerical 
techniques  were  employed  to  damp  or  filter  any  field.  The 
grid  refinement  obtained  when  using  a  rectangular 
subdivision  was  more  attractive  than  with  a  triangular 
subdivision.  It  allowed  refinement  to  be  easily  obtained 
with  Increased  accuracy  from  the  higher  order  polynomial  and 
decreased  computational  effort  due  to  reduced  operation 
counts.  Additionally,  It  was  shown  that  the  geometry  of  the 
triangular  elements  greatly  influences  the  problem. 
Furthermore  unless  specific  constraints  are  enforced,  the 
boundary  conditions  would  not  be  satisfied. 

Second,  the  ability  of  the  GFEM  to  resolve  atmospheric 
phenomena  occurring  near  the  scale  of  the  smallest  grid- 
length  was  shown.  A  small-scale  steady  state  solution  was 
cast  into  a  source  required  for  that  solution.  Integrations 
showed  that  the  GFEM  model  resolved  the  small-scale  forcing 
admirably  and  far  exceeded  the  performance  of  an  equivalent 
finite  difference  model.  Even  when  the  resolution  of  the 


H  finite  difference  model  was  doubled  and  tripled.  It  could 

not  match  the  performance  of  the  low-resolution  variable 
grid  GFEM  model.  The  responsiveness  of  the  GFEM  model  near 
the  smallest  resolvable  grid  length  clearly  demonstrates 
that  the  variable  grid  concept  Is  worth  the  computational 
effort. 

Finally,  the  concept  of  moving  the  variable  grid  was 
demonstrated.  The  ability  to  resolve  a  small-scale 
phenomenon  is  a  tremendous  asset  to  the  GFEM.  The  ability 
to  move  the  grid  with  the  atmospheric  phenomena  further 
enhances  Its  usefulness. 

Throughout  this  dissertation  the  GFEM  model  has 

I 

demonstrated  desirable  characteristics.  It  has  conservative 
properties.  It  propagates  atmospheric  waves  better  than  an 
equivalent  finite  difference  model.  It  allows  variable- 
resolution  grids  and  responds  better  than  an  equivalent 
finite  difference  model  near  the  smallest  gridlength. 

Moving  grids  can  be  achieved  with  no  apparent  noise 
generation.  It  can  utilize  direct  solvers  and  is  a  natural 
choice  for  vectori zati on  on  large  computers.  It  truly  is  a 
viable  option  for  simulation  of  atmospheric  flows. 

The  success  of  this  dissertation  mandates  further 
research.  The  next  logical  step  would  be  a  baroclinlc 
model.  Completion  of  a  baroclinlc  model  should  allow  its 
application  as  a  regional  model.  Additionally,  different 
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basis  functions  could  be  employed  for  the  vertical 
discretization,  such  as  solutions  to  the  vertical  structure 
equation.  For  the  Moving  variable-resolution  grid,  further 
research  Into  two-dl «ens1 onal  Movement  is  required.  The 
proof  of  concept  In  this  dissertation  of  one-dimensional 
movement  can  logically  be  extended  to  two  dimensions.  Upon 
completion  of  two-dimensional  movement,  an  operational 
forecast  capability  for  tropical  cyclones  should  be  Investi¬ 
gated.  The  superior  small-scale  response  of  the  GFEM 
indicates  potential  Increase  In  skill  for  both  regional  and 
tropical  cyclone/hurricane  prediction. 


APPENDIX  A 
STABILITY  ANALYSIS 


A  one-dl mensl onal  linear  analysis  of  the  forecast 
equations  Mill  be  presented.  The  analysis  was  obtained 
through  personal  communication  with  Dr.  Andrew  Stanlforth 
(Recherche  en  Prevision  Numerique).  The  primitive  form  of 
the  forecast  equations  and  a  semi -1 mpl 1  cl t  time  scheme  will 
be  analyzed,  because  the  results  are  Identical  for  the 
vortlclty/di vergence  form.  The  one-dimensional  equations 


with  a  mean  flow  U  are 


1H.+  11...U— +fv 
at  ax  sx 


OH  -  fu 
w  ax 


11  +  ,  .  0  11 

at  ^x  ax 


(A-l ) 


( A-2 ) 


(A-3) 


Time  derivatives  are  evaluated  with  a  centered  time  dif¬ 
ferencing  and  the  other  terms  on  the  left-hand  side  are 
averaged  between  time  levels  (t  +  at)  and  (t  -  at). 
Equations  (A-l),  (A-2)  and  (A-3)  become 


u(x. t+at)  -  u(x,t-at)  +  1 


,x+Ax,t+at)  -  a(x-ax.t+atl 
2ax 


[x+ax.t-at] 


[x-ax.t-at! 


■]  ■  -  u  C 


u(x+ax.t)  -  u(x-ax.t) 


+  f  v(x,t) 


(A-4) 


4  * 

1 


■Pi 


Assuming  that  a  function,  F,  can  be  expressed  as 
F(x,t)  -  F'e1(|tx 


then 


F (x A t+A t)^F(x, t-At),  a  IF*  s^n(wAt)e 


i(kx  +  wt) 


At 


F(«,t+at)  tF<x»*rAE)  -  F'  cos(wAt)e 


1 (kx  +  wt) 


Equations  (A-4),  (A-5),  and  (A-6)  become,  after 
substitution , 


{ A-7) 


(A-8) 


(A-9) 


1  su 1  , 
At 

fv*  +  1kc$'  *  0 

( A— 1 0 ) 

fu'  + 

hi'  ,  o 

At 

(A-m 

1kc«u' 

+  In'  .  o 

At 

( A- 12) 
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where  s  ■  sin  wAt  +  k'uAt,  c  ■  cos  wAt  and  k'  *  sin 
The  determinant  of  the  system  of  equations  (A-10),  (A-11) 
and  (A-12)  when  set  equal  to  zero  yields 

slr^-tf^kV.n.o 

The  roots  are  s*0  and 

s2  -  IfAt)2  +  c2CkAt)2  *  (A"14 

The  roots  of  this  second  order  equation  when  requiring  w  to 
be  real  yield  the  stability  criterion 


At 


|-^|  +  f 
AX' 


(A-15) 
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APPENDIX  B 
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i 


i 

4 
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NUMERICAL  QUADRATURE 

Evaluation  of  the  Galerkin  Integrals  requires  a  fast, 
efficient  method  .  Two  different  elements,  triangles  and 
rectangles,  have  been  employed  and  the  quadrature  schemes 
for  evaluation  are  presented  here.  In  the  first  case, 
triangles,  area  coordinates  are  used.  For  the  rectangles, 
integration  formulas  are  based  on  an  orthogonal  axis  trans 
formation.  In  both  cases,  the  integrals  to  be  evaluated 
contain  either  products  of  basis  functions,  products  of 
derivatives  of  basis  functions,  or  a  mixture  of  both. 

Zienkiewicz  (1971)  developed  the  relationship  between 
the  Cartesian  coordinates  of  a  triangle  and  the  area 
coordinates  as 


X  *  L1X1  +  L2X2  +  L3X3 
Y  =  L,Yl  +  L2Y2  ♦  L3Y3 


(B-l) 

(B-2 ) 

(B-3) 


where  (xI,y1)v  (x2,y2)  and  (x3,y3)  are  the  Cartesian  coordinates 
of  the  triangle's  vertices  and  L^,  l2  and  L3  are  the  area 
coordinates.  Here 

L1  »  A,/A  l2  «  A2/A  L3  »  A3/A 
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where  A^,  and  ^3  a1*®  shown  In  Figure  54  and  A  is  the  area  of 
the  entire  triangle.  He  further  shows  that 


//ha  4b  l3C  dxdy 


a.'b.'c! 


The  L's  are  the  basis  and  test  functions  used  during  the 
approximating  process.  Figure  55  defines  some  necessary 
parameters  for  evaluation  of  the  derivatives. 

Equations  B-l,  B-2  and  B-3  can  be  written  in  matrix 


form  as 


2A  b1  a-jl  rll 


L2  =  2A  2A  b2 
L~  2A  b, 


a3JLY- 


(B-5) 


Differentiation  of  Eq.  (B-5)  shows  ♦'hat 


_La  l  bi  a 
3x  ~  .=1  2A  3Lt 


( B  —  6 ) 


_L-  I  !i_i_ 

ay  ii1  2A  31. 


( B  -  7 ) 


However,  the  derivative  of  the  basis  function  with  respect 
to  the  area  coordinate  is  non-zero  only  where  the  basis 
function  and  area  coordinate  coincide.  In  fact,  the  non¬ 
zero  value  is  one.  Therefore,  using  Vj  as  a  basis  function 


3V  i  A 

3x“  2A 


(B-8) 


All  integrals  can  be  evaluated  by  combination  of  these 
quadrature  formulas  for  triangles. 

An  orthogonal  axis  transformation  will  allow  similar 
quadrature  formulas  for  rectangles.  The  rectangle  shown  in 
Figure  56  can  be  transformed  by 


The  values  of  ?  and  n  at  each  corner  are  shown  in 

parentheses.  The  basis  function,  ,  becomes 

(l+Ci?)(l+nin) 
vi  4 


1-1.11 


(XoJcl 


1-1.-1 1 


11.-1) 


Fig.  56.  Orthogonal  axis  transformati on  for  rectangular 
integration  formulas. 

Derivatives  of  the  basis  function  become 


3 


By  substituting  and  Integrating,  the  4x4  matrix  with  the 
interaction  coefficients  can  be  determined  for  a  derivative 
or  straight  inner  product.  For  instance,  the  straight  inner 
product  1 s 

ctj  ■//vivjdxd^  ■ ab  1 1  Vj  d'd'> 

i  i 

■  J  (l+n1n)(l+njTi)dn 


and  the  mixed  derivative  is 
•3  V,  3V, 


I  =(f- 1 
ij  JJ  3x 


1  1 


3X 


dxdy 


-  f  f  -1 

*  /j  /j  3? 


3V,  3V 


ac 


dcdn 


i  i 

m  1 5tcjdi;  /,  ll4,'i")(,+v)d" 


_*  "iV  (B‘U) 

These  quadrature  rules  allow  the  evaluation  of  any  integrals 
when  using  rectangles. 
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